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We carry out extensive Monte Carlo simulations of the three-dimensional (3D) uniformly frus- 
trated XY model with uncorrelated randomly perturbed couplings, as a model for the equilibrium 
behavior of an extreme type-II superconductor with quenched uncorrelated random point vortex 
pinning, in the presence of a uniform applied magnetic field. We map out the resulting phase di- 
agram as a function of temperature T and pinning strength p for a fixed value of the vortex line 
density. At low p we find a sharp first order vortex lattice melting phase boundary separating a 
vortex lattice from a vortex liquid. As p increases, it appears that this first order transition smears 
out over a finite temperature interval, due to the effects of the random pinning, in agreement with 
several recent experiments. At large p we find a second order transition from vortex liquid to vortex 
glass. 
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I. INTRODUCTION 

In strong type-II superconductors, where the magnetic 
penetration length A is much larger than the bare coher- 
ence length £o> much of the macroscopic behavior can be 
modeled in terms of the interacting vortex lines that are 
introduced into the system by the application of an exter- 
nal magnetic field. 1 In the high temperature copper-oxide 
superconductors such as YBCO and BSCCO, material 
parameters such as high anisotropy, as well as the high 
transition temperatures, cause thermal fluctuations to 
play an important role. Such materials are also believed, 
even in the purest of single crystal samples, to contain 
quenched, intrinsic, uncorrelated, random point impuri- 
ties that can pin vortex lines. The resulting phase dia- 
gram as a function of temperature and applied magnetic 
field is generally believed to result from such a combina- 
tion of thermal fluctuations and random point pinning. 2 
It was later argued that similar fluctuation effects can be 
observed even in single crystal samples of more familiar 
strong type-II low temperature superconductor^ 3 -* 4 1 5 1 6 1 7 1 
such as Nb and NbSe2, albeit over a much more restricted 
region of the phase diagram. It has been argued that a 
universal phase diagram might apply to all strong type-II 
superconductors P 

Considerable experimental effort has been devoted 
to the determination of this vortex line p hase 
diagram.^' 11 ' 12 ' 13 ' 14 ' 15 ' 16 '^^ is 
now generally accepted that upon increasing temperature 
at low magnetic fields, a sharp first order melting transi- 
tion exists^ ^ 1 * 12 * 13 * 14 ! from an elasticall y disto rted vortex 
lattice, known as the "Bragg glass, ;, HHl 29 l 3Q l to a disor- 
dered vortex liquid. However many other aspects of the 
vortex phase diagram remain in dispute. At low temper- 
atures, increasing the magnetic field leads to a first order 
trans ition from the Bragg glass to a disordered vortex 
statePUSl it was originally believed that a special criti- 



cal point separated the low field thermally induced vortex 
lattice melting fr om the high field disorder ( random pin) 
induced meltingPEMI] However later wor M 15 ' 16 ' 17 ' ^ 
argued for a single unified first order phase boundary 
for the vortex lattice, continuing smoothly down to low 
temperatures. At low temperatures and high magnetic 
fields, where the vortex lines are spatial ly dis ordered, a 
"vortex glass" phase has been proposed! 31 * 32 * It remains 
in question whether this vortex glass is a truly distinct 
thermody namic state with true superconducting phase 
coherence, 1 22 1 23 1 separated from the vortex liquid by a con- 
tinuous second order phase transition, or whether there 
is just a crossover to a highly viscous vortex liquid upon 
decreasing temperature! 24 1 25 1 It has also been proposed 
that within the vortex liquid there is a sharp first or- 
der transition line that splits off from the melting line at 
higher magnetic fields, and terminates at a critical end 
point. This transition is claimed to separate regions with 
lesser vs greater spatial vortex correlations, th e mo re cor- 
related region being called the "vortex slush" J^ 2 ^ Other 
works have argued that the disordered vortex state at 
high magnetic fields continues to exist as a thin sliver 
of "multidomain glass" all along the vortex melting line, 
even at low magnetic fields 

Theoretically, the vortex line phase diagram has been 
studied within effective elastic theories, often using the 
Lindemann criterion to e stimate the lo cation of melt- 
ing a nd othe r transitions! 33 | 34 | 35 | 36 | 3T | 38 1 Several of these 
work d 36 ! 37 ! 38 ! reported the possibility of a critical end 
point to the vortex lattice melting line at high magnetic 
fields, as well as a vortex slush phase. However later 
analytical calculations, based on the Ginzburg-Landau 
model with pinningj 3 -^ and on a vortex line model in a 
weak impurity background with both elastic and plas- 
tic excitations, 40 found only a unified first order melting 
transition for the Bragg glass, extending to low temper- 
atures with no critical points intervening, and no vortex 
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slush. 

The difficulty of performing reliable analytical 
calculations has thus led to numero us numerical 

in which one 

hopes to map out various aspects of the vortex line phase 
diagram within a well defined simple model. In this pa- 
per we report on one such investigation, using the three- 
dimensional (3D) uniformly frustrated XY model with 
uncorrelated quenched disorder in the couplings, as a 
model for a strong type-II superconductor with random 
point pinning. A virtue of the XY model is that it de- 
scribes realistic vortex line interactions with no restric- 
tions on vortex line excitations; both overhangs in the 
field induced vortex lines, and closed thermally excited 
vortex loops are included. The details of vortex line cores 
are handled by the short length cutoff of the numerical 
grid, and hence issues concerning vortex line cutting are 
treated with a minimum of ad hoc assumptions. Our 
model holds in the limit of infinite penetration length, 
A — » oo, and we comment later on the implications of 
this approximation. 



We build upon our earlier worlP^^ to present here 
the equilibrium phase diagram as a function of temper- 
ature T and disorder strength p, for a fixed density of 
vortex lines / = B^/fa = 1/5 (</>o is the flux quantum). 
Keeping / fixed avoids effects that would be due to vary- 
ing commensurability of the vortex lines with respect to 
the underlying numerical grid of the model. However, 
since increasing / in a continuum system is generally 
believed to increase the effective pinning strength,^ as 
more lines get forced into the same pinning volume, vary- 
ing p at fixed / should provide qualitatively similar in- 
formation as the more physical situation of varying / 
at fixed p. Through careful, well equilibrated, simula- 
tions comparing systems of different size, and averaging 
over different realizations of the random pinning, we con- 
sider the limits of both weak and strong disorder. Our 
results are summarized in Fig. [I] By using a denser vor- 
tex line system than in our earlier work we are able in 
particular to explore the strong pinning limit, correct- 
ing our earlier preliminary conclusions^ concerning the 
presence of a vortex glass in this model. We also find, for 
the first time, a disorder induced smearing of the vortex 
lattice melting transition T m (p) at intermediate disorder 
strengths, that w e be lieve is in good agreement with re- 
cent experiments F^U 

The remainder of this paper is organized as follows. In 
Sec. II we define our model and the quantities we use to 
map out the phase diagram. In Sec. Ill we present our 
results for vortex lattice melting in the limit of weak pin- 
ning strengths. In Sec. IV we present our results for the 
vortex glass at strong pinning strengths. In Sec. V we 
discuss the smearing of the vortex lattice melting tran- 
sition at intermediate pinning strengths. In Sec. VI we 
summarize our results and discuss their relation to other 
recent numerical works. 
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FIG. 1: (Color online) a) Phase diagram of the point disor- 
dered 3D XY model of Eq. Q as a function of temperature T 
and disorder strength p for fixed vortex line density / = 1/5 
and coupling ratio J z = J±/40. The solid blue circles (•) and 
solid line at low p indicate a sharp first order vortex lattice 
melting transition T m (p) , as determined from a system of size 
L = 20, L z = 6 and averaged over 8 independent realizations 
of the quenched randomness. As p increases, this melting 
transition is observed to smear out over a finite temperature 
interval, T m 2 < T < T m i, as determined from a system of size 
L = 20, L z = 12 and averaged over 4 — 8 independent real- 
izations of the quenched randomness; T m i is indicated by the 
open red squares (□) and dashed line, while T m 2 is indicated 
by the open purple circles (o) and dashed line. Error bars 
represent the standard deviation of values obtained over the 
independent realizations. At large p we find a sharp second 
order transition T g (p) from vortex liquid to vortex glass, as 
determined by a detailed finite size scaling analysis of system 
sizes L = 10 - 20, L z = (3/5)L, averaged over 200 - 600 
independent realizations of the quenched randomness; T s (p) 
is indicated by the solid black squares (■) and solid line. 
We are unable to equilibrate our system at the low tempera- 
tures (denoted by "???") where the vortex melting and vortex 
glass transitions approach one another, b) Expanded view of 
the phase diagram showing the smearing of the vortex lattice 
melting transition as p increases. 
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II. MODEL 

The model we simulate is the uniformly frustrated 3D 
XY model, given by the HamiltonianJ^ 

= - J v* cos (#(r*) - 0(v % + A) - Aip) , (1) 



which serves as a discretized approximation to the 
Ginzburg-Landau free energy functional in the London 
approximation of fixed wavefunction amplitude, and cap- 
tures the kinetic energy of the flowing supercurrents in 
the system. Here Ofa) is the thermally fluctuating phase 
angle of the superconducting wavefunction on site of 
a cubic L x L x L z grid of sites with bonds in directions 
pi = x,y : z. J iyb is the superconducting coupling constant 
on bond (i/i) of the grid, A ifl is proportional to the mag- 
netic vector potential integrated across the bond, and the 
argument of the cosine is thus the gauge invariant phase 
angle difference across the bond. 

The are determined such that the circulation of 
A ifI around any plaquette a of the grid is fixed and equal 
to 27r/ a with f a equal to the fraction of applied magnetic 
flux quanta through that plaquette. To model a uniform 
applied magnetic field in the z direction, we use a uni- 
form value f a = f for plaquettes oriented with normal in 
the z direction, and f a = otherwise. The presence of 
the applied magnetic field induces vortex lines into the 
phase angles 0(ri) such that the average vortex line den- 
sity in direction fi is n M = f5^ z . Holding the fixed 
corresponds to the approximation of an infinite magnetic 
penetration length A, which one may expect to be reason- 
able in the limit that A is much larger than the average 
inter-vortex spacing. Experiments on very pure YBCO 
sampled indicate that many features of behavior, up to 
surprisingly large magnetic fields, are well described by 
this large A — > oo XY model limit. We return to comment 



on this approximation in section IV 



To model vortex pinning due to quenched point ran- 
domness we use couplings 46 

Jifi = H = x,y (2) 

Jiz = Jz constant, 

where the e iyb are uncorrelated, uniformly distributed, 
random variables with 



M = 0, (e? ) = 1 . 



(3) 



The disorder strength is thus controlled by the parameter 
p in Eq. (J2|. We will vary p from small to large values 
in order to systematically investigate the differences be- 
tween weak and strong pinning. 

In the following we will use (...) to denote the equilib- 
rium average for a particular realization of the random 
pinning {e^}. We will use [. . .] to denote the average 
over several independent realizations of the {e^}. 

In this work we use parameters, 



/ = 1/5 , 
Jz = J±/40 . 



(4) 
(5) 



We use J z <C J± to enhance vortex line fluctuations along 
the z direction, thus allowing us to use systems with 
smaller L z . The relatively dense value of / was similarly 
chosen so as to have many vortex lines contained within 
systems of modest size, so as to permit finite size scaling 
analyses (particularly for the vortex glass) and to average 
over many independent realizations of the random disor- 
der {e^}. While this high density leads to artificial com- 
mensurability effects (for example a square rather than 
a triangular ground state vortex lattice), one may still 
hope that many features of our system with respect to 
vortex lattice melting and the vortex glass transition will 
remain qualitatively similar to results on real physical 
materials. 

To probe the behavior of the system we consider the 
average energy per site, which is the thermodynamic con- 
jugate variable to the temperature T, 



E 



1 d(/3T) 



L 2 L Z 



L 2 L 



dp 



(6) 



(cos(0(r^-0(r; + /i)-^)) 



where (3 = 1/T and T is the total free energy. We can 
similarly define a variable Q which is the thermodynamic 
conjugate to the disorder strength p, 

Q s ik% m 

= ~JJl~ H J±e ill (cos(0(r i )-0(r i + ji)-A ill )) . 

The vorticity in the system is determined by consider- 
ing the circulation of the gauge invariant phase difference 
around each plaquette. For a plaquette a at position r a 
with normal in direction /i, the vorticity n M (r Q! ) piercing 
the plaquette is determined by, 



[0(n) - #(r, + fi) - Ai^Zl = 2^[n M (r a ) - fS^] , 



(8) 

where the sum goes counterclockwise around the bonds 
of the plaquette, and the bracket on the left hand side 
indicates that the gauge invariant phase angle difference 
is to be computed so as to lie within the interval (— 7r, tt]. 

To look for vortex lattice ordering we compute the vor- 
tex structure function 5(1c_l) which measures correlations 
between vortices within the same xy plane, 



5(k ± ) 



fL 2 L z 



e ik -( r - - r '-) (n z (rx , z)n z (r' ± , z) ) 

(9) 

where = (x, y) denotes the coordinates in the xy plane 
and similarly = (k Xl k y ). The presence of a vortex 
lattice will be indicated by the appearance of sharp peaks 
in £(k_L) at reciprocal lattice vectors K. 

To look for a possible vortex glass phase, in which vor- 
tex lines are frozen into a disordered configuration and 



4 



so 5(kj_) displays no signal of the ordering, we consider 
the helicity modulus.^ To simulate the Hamiltonian of 
Eq. on a finite size grid, we use fluctuating twist 
boundary conditions, defined by, 



0(r< + Lf.fi) - 0(n) = A p 



(10) 



where A M , the total phase angle twist across the system 
in direction /}, is taken as a thermally fluctuating degree 
of freedom. We then compute the histogram P(A M ) of 
values that A M takes during the course of the simula- 
tion (averaging over the twists in the transverse direc- 
tions) and define the free energy variation with twist by, 
^(A^) = — TlnP(A /i ) + constant. The helicity modulus 
T M in direction [i is then defined in terms of the curvature 
of ^(A^) at its minimum A 



T = 

1 (J- — 



/xOj 



L V L 



dAl 



(11) 



where /i, v, a are a permutation of x, y, z. Note, unlike 
pure systems, for random systems it is not generally true 
that A^o = 0. When ^(A^) varies with A M , and so 
the system is sensitive to the boundary conditions, we 
have T M > and the system possesses superconducting 
phase coherence. When ^(A^) is flat and independent 
of A M , T M = and superconducting phase coherence 
is lost. The vanishing of T M thus is a signature of the 
superconducting transition. 

At a second order phase transition T c , the tempera- 
ture and system size dependence of T M is expected to 
obey a critical scaling law. Since the applied magnetic 
field singles out a special direction, there is the possibil- 
ity that this scaling may be anisotropic. In such a case, 
the expected scaling law for T M (T, L,L Z ) is, 

%^T M (T, L, L z ) = u^tl}l\ L x /L<) , (12) 

where t = T — T c , u is the scaling function, v is the cor- 
relation length critical exponent, and £ is the anisotropy 
critical exponent. Should the scaling turn out to be 
isotropic, then £ = 1, and for systems with a fixed as- 
pect ratio L z = jL the scaling law reduces to, 



LT ll {T,L)=u tl {tL 1 l v ) 



(13) 



In this case, exactly at the transition temperature T c , one 
has t = and so LT M is independent of system size L. 

Henceforth, we will measure temperature and energies 
in units where J± = 1. Length will be measured in units 
where the grid spacing is unity. 



III. VORTEX LATTICE MELTING AT WEAK 
PINNING 

In this section we consider the first order vortex lat- 
tice melting transition at weak pinning strength. Our 
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FIG. 2: a) Real space configuration of vortex lines (denoted 
by "+" ) piercing the xy plane in the ground state of the pure 
(p — 0) model, for a vortex line density / = 1/5. The two 
possible orientations of this ground state vortex lattice with 
respect to the underlying numerical grid are shown, b) In- 
tensity plot of vortex structure function <S , (kj_) for each of 
the two ground state vortex lattice orientations shown in (a); 

= is at the center of the plot. Peaks at the non-zero 
reciprocal lattice vectors are labeled Ki and K 2 respectively 
for the two orientations. 



methods for identifying the melting transition and es- 
tablishing that it is indeed a first order phase transition 
are the same as we have used in our earlier workPD on the 
more dilute / = 1/20 system. 

For a pure system with disorder strength p = the vor- 
tex line lattice, of vortex density / = 1/5, will order into 
a ground state containing a square vortex lattice with 
lattice constant y/E. There are two possible orientations 
of this square lattice with respect to the underlying grid, 
related to each other by a reflection through the x axis, 
as shown in Fig. [2^i. In Fig. we show the vortex struc- 
ture function £(kx) for each of these two ground state 
orientations. We find sharp Bragg peaks at reciprocal 
lattice vectors K. For the vortex lattice of Fig. [2^i there 
are only four non-zero reciprocal lattice vectors, related 
to one another by tt/2 rotations of /c-space. For the two 
possible ground state orientations, we label these two dis- 
joint sets of reciprocal lattice vectors by Ki and K2 as 
shown in Fig. [2]o. 

At finite but small disorder strength, p < p c c± 0.22, 
we continue to find at low temperatures a vortex lat- 
tice with the same symmetry as that of Fig. |2^i. The 
existence of the two possible orientations for this state 
motivates our definition of the following order parame- 
ter for the vortex lattice to liquid melting transition. If 
we denote by £(Ki) and S(K.2) the value of the vortex 
structure function averaged over the four non-zero recip- 
rocal lattice vectors of each orientation respectively, then 
the difference, AS = 5(Ki) — 5(K2), will signal the vor- 
tex lattice melting transition: below melting, the system 
has ordered into one of the two possible vortex lattice 
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orientations, and so S(Ki) is large and S(K 2 ) is small, 
or vice versa, giving a large value of |AS|; above melting, 
the system is in a liquid state with the same symmetry 
as the underlying grid, so S(Ki) = S^K^) by reflection 
symmetry and AS = 0. 

Anticipating a first order vortex lattice melting tran- 
sition, we slowly cool down a 20 x 20 x 6 size system 
from high temperature using ordinary Metropolis Monte 
Carlo until we reach a temperature at which we observe 
the system to switch back and forth between large and 
small values of AS during the course of the simulation. 
As an example of this, we show in Fig. Km a plot of AS/ So 
vs. simulation time (So = S(0) = fL 2 ) at a tempera- 
ture close to the melting transition T m (p) for a moderate 
value of disorder strength p = 0.12. In Fig. ^Sjp we show 
an intensity plot of lnS(kj_), averaged over only those 
configurations in Fig. [3^ which have AS /So > 0.1. We 
see sharp Bragg peaks of high intensity at the recipro- 
cal lattice vectors Ki indicating that this is the vortex 
lattice state. In Fig. ^ we show show an intensity plot 
of lnS(k^), averaged now over only those configurations 
in Fig. [3^ which have AS /So < 0.1. We see broad dif- 
fuse peaks of equal low intensity at both Ki and K2, 
indicating that this is the vortex liquid state. 

In Fig. [4^ we plot the histogram P(AS) of values of 
AS found in the data of Fig. [3^l. We see two well sepa- 
rated peaks centered at AS /So — and AS /So — 0.55, 
representing the vortex liquid and vortex lattice states 
respectively. We define the melting transition temper- 
ature T m (p) to be the temperature at which the area 
under these two peaks is equal. In this manner, vary- 
ing p and averaging over 8 independent realizations of 
the random disorder, we plot the vortex lattice melting 
line T m (p) for systems of size 20 x 20 x 6 as the blue solid 
curve m Fig. 1 As p increases, r m (p) decreases while the 
slope \dT m /dp\ rapidly increases. As p increases towards 
p c ~ 0.22, with correspondingly low melting T m , a failure 
to achieve proper equilibration of the system prevents us 
from continuing to trace out the melting curve to lower 
temperatures. 

Next we demonstrate that, within our model system, 
melting remains a first order transition along the melting 
curve for as far as we can map it out. There is no sign 
of it ending at an "upper critical point" as has been of- 
ten suggested by experimental works Q^EIl Choosing the 
minimum in the P(AS) histogram as the dividing point, 
we assign each configuration as a vortex lattice or vortex 
liquid according to the value of AS for that configura- 
tion. Having divided configurations into distinct lattice 
and liquid states, we can then construct the histograms 
of energy, P(E), and of the disorder conjugate variable, 
P(Q), for each state respectively. In Figs. ^jp,c we show 
such histograms for disorder strength p = 0.12, corre- 
sponding to the data of Figs. [3]and[4^L. We see for both 
E and Q well separated histograms for lattice and for 
liquid states. Using these histograms we then compute 
the average E and Q separately for the vortex lattice 
and vortex liquid, and then compute the discontinuities 
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FIG. 3: a) Plot of AS = S(Ki) - S(K 2 ), normalized by 
So = S(k^ = 0), vs Monte Carlo simulation time, close to 
the first order melting transition temperature T m (p) c± 0.228 
for a particular realization of the quenched randomness at 
disorder strength p = 0.12. The unit of time on the plot rep- 
resents 10 6 Monte Carlo sweeps through the entire system of 
grid size L = 20, L z = 6, and each data point represents an 
average over 65536 consecutive sweeps. We see the system 
making many discrete transitions between vortex liquid, rep- 
resented by small values of AS /So ^0.1, and vortex lattice, 
represented by large values of AS/ So > 0.1. b) Logarithmic 
intensity plot of vortex structure function lnS(k^) averaged 
over only the vortex lattice states of (a). A sharp Bragg peak 
at Ki is seen, c) Logarithmic intensity plot of lnS(kj_) aver- 
aged over only the vortex liquid states of (a). Diffuse peaks 
of equal height are seen at both Ki and K 2 . 

in these quantities at melting, 

AE = (illiquid - (^lattice (14) 
AQ = (Q) lattice — (Q)liquid 

In Fig. |5|i we plot (open symbols) the resulting [AE] 
and [AS] , averaged over 8 independent realizations of the 
random disorder, vs disorder strength p, for the system 
of size 20 x 20 x 6. As p increases we see that [AE] de- 
creases and appears to vanish. A vanishing [AE] implies 
a vanishing entropy jump and hence a vanishing of the 
delta-function specific heat singularity usually associated 
with the first order melting transition. Such a vanish- 
ing of the specific heat singularity, observed experimen- 
tally on increasing applied magnetic field at fixed disorder 
strength, has been used as evidence for a weakening of 
the first order melting transition and its termination at a 
second order upper critical point However, if the first 
order phase transition is indeed to vanish at a second or- 
der critical point, it is necessary that the discontinuities 
in all thermodynamic first derivatives of the free energy 
vanish as the critical point is approached. In Fig. [5^i, 
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FIG. 4: (Color online) a) Histogram P(AS) of values of 
AS /So obtained at the melting temperature T m (p) for a par- 
ticular realization of the quenched randomness at disorder 
strength p = 0.12, system size L = 20, L z — 6, from the data 
of Fig. [3^ . The minimum at AS /So ~ 0.1 separates config- 
urations corresponding to the vortex liquid, represented by 
the peak centered at AS /So = 0, from the configurations 
corresponding to the vortex lattice, represented by the peak 
centered near AS /So ~ 0.55. b) Histograms P(E) of energy 
per grid site E, and c) histograms P(Q) of the disorder conju- 
gate variable Q, for configurations in the vortex liquid phase 
(AS/ So < 0.1) and the vortex lattice phase (AS/S > 0.1) 
of (a). Histograms indicate discontinuous jumps AE and AQ 
between vortex liquid and lattice at the melting transition. 



however, we see that the discontinuity [AQ] increases as 
p increases, and does not vanish as it must if the first or- 
der line is to end in a critical point. We thus find that our 
melting transition remains strongly first order for all dis- 
order strengths p. The vanishing of [AE] upon increasing 
p merely reflects the increasing slope of the melting curve 
\dT m /dp\ in accordance with the Clausisus-Clapeyron re- 
lation, as can be seen as follows. Since the free energies 
of lattice and liquid must be equal at T m , 

AF[T m (p),p] = lattice [T m (p),p]- liquid [T m (p),p] = , 

(15) 

we have, 



dAF dAF dAF dT„ 



dp 



dp 



dT dp 



Since 
dAT 



dp 



L L Z AQ and 



dAf 
dT 



L Z L Z 



AE 



(16) 



(17) 



substituting into Eq. (|T6|) then gives the Clausius- 

(18) 



Clapeyron relation for our system, 
dT m T m AQ 



Fitting our data (blue solid circles in Fig.JT]) for T m (p) to 
a quadratic polynomial in p 2 (solid blue line in Fig. [I]), 
we used the fitted polynomial to determine the slope 
dT m /dp, and in Fig. [5b we plot \dT m /dp\ and the dis- 
order averaged [T m AQ/A^] vs disorder strength p. The 
disorder average is over 8 independent realizations of the 
quenched randomness. We find excellent agreement with 
Eq. (18), thus verifying that our results are indeed very 



dp 



AE 



well equilibrated. 

Experiment evidence for the absence of an upper crit- 
ical point, in agreement with our results, has been ob- 
tained in BSCCO by Avraham et al. 17 For an exper- 
imental system, in which disorder strength is constant 
and the applied magnetic field H is varied, the magne- 
tization density M = (l/V)dJ r /dH becomes the analog 
of our parameter Q. While initial measurements of the 
jump AM at melting appeared to show AM vanishing 
as H increased, suggesting an upper critical point, subse- 
quent measurements using a additional small oscillating 
field to "tickle" the vortex lines to help avoid trapping 
in metastable local energy minima, showed a finite AM 
continuing along the melting curve past the presumed up- 
per critical point and down to even lower temperatures. 
Their conclusion was that the presumed upper critical 
point in BCSSO was an artifact of poor equilibration, 
and that a unified first order transition line continued 
between thermally driven melting at low H, and disor- 
der driven melting at large r H. Similar conclusions had 
been drawn earlier by others^MJ based on measurements 
of the Josephson plasma frequency. 

Finally, we consider the finite size dependence of our 
results to see that the values of [AE] and [AQ] which we 
have found for the 20 x 20 x 6 size system do not appear to 
be decreasing (or perhaps vanishing) as the system size 
increases. The need to keep the transverse system length 
L a multiple of 5, so as to remain commensurate with the 
vortex lattice periodicity, and the difficulty of equilibrat- 
ing hops between lattice and liquid states as the system 
size, and hence to the total free energy barrier between 
theses state, increases, limits greatly the range of system 
sizes that we can consider. In Table [J we show our results 
for the three system sizes 20 x 20 x 6, 30 x 30 x 6 and 
20 x 20 x 12, for the specific disorder strength p = 0.12. 
We see that while T m decreases slightly as the system 
size increases, [AE] remains remain roughly independent 
of size while [AQ] shows a slight increase. We also com- 
pute the spread in melting temperatures AT m that we 
find as we consider different independent realizations of 
the quenched random disorder {e^}. If the system is self 
averaging over the quenched disorder, we would expect 
that AT m oc l/y/V, with V = L 2 L Z the system volume. 
In Table [i] we therefore also give the value for AT m \^V 
for the three system sizes. Considering the relatively few 
(i.e. 8) disorder realizations we have considered for the 
two larger sizes, and hence the corresponding large po- 
tential error in our estimate of AT m , we find our results 
consistent with the expectation of self averaging. 
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FIG. 5: (Color online) a) open symbols: Disorder averaged 
discontinuities in the energy per site [AE] and disorder conju- 
gate variable [AQ] at the melting transition T m (p) vs disorder 
strength p, for systems of size L = 20, L z = 6; solid symbols: 
average total changes [AE] and [AQ] in ordering from liquid 
to solid for systems of size L = 20, L z = 12 (see discussion 
Sec. \v\ . While [AE] decreases with increasing p, [AQ] re- 
mains finite and grows, indicating that the melting transition 
remains strongly first order, b) Comparison between the slope 
of the melting phase boundary dT m /dp and the disorder aver- 
aged value of \T m AQ/AE] vs disorder strength p, as a test of 
the Clausius— Clapeyron relation, using data for system size 
L = 20, L z = 6. Good agreement indicates that our results 
are well equilibrated. Error bars represent the estimated sta- 
tistical error sampling over independent realizations of the 
random disorder. 



TABLE I: Disorder averaged melting temperature [T m ] and 
discontinuities in the energy per site [AE] and disorder con- 
jugate variable [AQ] for disorder strength p = 0.12, for dif- 
ferent system sizes L 2 x L z . The errors represent the esti- 
mated statistical error sampling over the number of indepen- 
dent realizations of the random disorder specified in the last 
column. AT m is the standard deviation of melting temper- 
atures computed over the independent random realizations, 
and one expects AT m oc 1/y/V, with V — L 2 L Z the volume 
of the system. Results for the 20 2 x 12 system are computed 
as described in Sec. IVl 



L 2 x L z 




[AE] 


[AQ] 


AT m VV 
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20 2 x 6 


0.2304 
± 0.0009 


0.028 
± 0.0005 


0.057 
± 0.002 


0.183 


20 


30 2 x 6 


0.227 
± 0.001 


0.026 
± 0.002 


0.061 
± 0.004 


0.190 
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20 2 x 12 


0.225 
± 0.001 


0.0279 
± 0.0004 


0.067 
± 0.001 


0.184 
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IV. VORTEX GLASS AT STRONG PINNING 

We now consider the strong pinning limit, p > p c ~ 
0.22, and the possible existence of a vortex glass phase. 
In our earlier work 46 on the more dilute / = 1/20 sys- 
tem, we presented preliminary evidence for the absence of 
a vortex glass phase within the XY model. However later 
work by OlssorP^ established that the model does indeed 



have a second order vortex glass transition at strong dis- 
order strength. Moreover he found that critical scaling 
is isotropic. Here we follow the analysis of Olsson, look- 
ing at the disorder averaged helicity modulus in the xy 
plane, [T_J, defined in Eq. ( pTj ). 

In Fig. [6|i we plot [XTJ vs T for three different dis- 
order strengths, p = 0.3,0.4 and 0.55, using for each 
case three different system sizes from L = 10 to 25, 
L z = (3/5)L, as indicated in the figure. Results are 
averaged over 200 — 600 independent realizations of the 
random disorder, depending on system size. As discussed 
following Eq. ( fl3j ) , for each value of p, the common in- 
tersection of the curves for different L locates the vortex 
glass transition, T g (p), and thus allows us to map out the 
vortex glass transition line, plotted as the black squares 
and black solid line in Fig.[l] We see that T g (p) increases 
for increasing p. 

In Fig. ^jp we replot our data vs the scaled temperature 
tL x l v ', where the critical exponent v has been chosen for 
each p so as to give the best data collapse to a common 
scaling curve for the different system sizes L, in accor- 
dance with the scaling equation of Eq. (13). For the two 
largest values p = 0.4 and 0.55 we find v ~ 1.5, in agree- 
ment with the more precise calculations of Olsson. For 
the smallest p = 0.3, we find a smaller v ~ 1.3. We be- 
lieve that for this last case, where our data (see Fig. [gJl) 
is the noisiest and as T g is the smallest we have the least 
data at T < T g , the fact that we are closest to the vor- 
tex lattice melting transition may mean that our system 
sizes are still too small and so we are in a cross-over re- 
gion rather than the true large L scaling limit. 

The existence of a vortex glass phase in the bond dis- 
ordered 3D uniformly frustrated XY model has also been 
claimed in simulations by Kawamuraj 48 * 49 * who in addi- 
tion to measuring a value v = 1.1 ±0.2 also measures the 
critical exponent r] = —0.5 ±0.1. These are also consis- 
tent with the values v = 1.3±0.2 and r] = — 0.4±0.1 found 
by Lidmar 5 -^ using an elastic model of interacting vortex 
lines that includes dislocations. These values for the crit- 
ical exponents lie in reasonable agreement with the values 
v GG = 1.39 ± 0.20 and rj GG = -0.47 ± 0.07 found for the 
simpler 3D gauge glass model! 5 - 5 ^ The gauge glass has the 
same Hamiltonian as Eq. ([I]), except that the bond cou- 
plings Ji^ are taken as uniform, and the randomness is 
put in the vector potential, with A ifI chosen randomly 
from a uniform distribution on (— 7r, tt]. In the gauge 
glass, the average magnetic field thus vanishes in all di- 
rections and so the model is intrinsically isotropic. This 
comparison of critical exponents suggests that the vortex 
glass and the gauge glass may be in the same universality 
class. If so, Kawamura notes 49 ' that the addition of mag- 
netic field fluctuations associated with a finite magnetic 
penetration length A ( "magnetic screening" ) in a real su- 
perconductor may serve to destabilize the vortex glass 
transition and replace it with a smooth crossover behav- 
ior, as has been numerically observed to happen in the 
gauge glass 56 (such an effect was observed in early simu- 
lation by Kawamura 4 ^ of the vortex glass with magnetic 
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V. SMEARED MELTING AT INTERMEDIATE 
PINNING 

We return now to the vortex lattice melting at p < 
p c ~ 0.22. Considering a system size 20 x 20 x 12, twice as 
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FIG. 6: (Color online) a) Transverse helicity modulus scaled 
by system length LT^ vs temperature T for different sys- 
tem sizes L, L z = (3/5) L, for strong disorder strengths 
p — 0.30, 0.40 and 0.55. The common intersections of the 
curves for different L locate the vortex glass transition tem- 
peratures T s (p). Results are averaged over 200 independent 
realizations of the random disorder for the largest system size, 
and between 300 — 600 realizations for smaller sizes, b) LT^ 
of (a) replotted vs scaled temperature tL 1 ^ , with t = T — T g . 
Collapse of the data for different L to a common scaling curve 
determines the correlation length critical exponent v ~ 1.5. 



screening, however the use of free boundary conditions 
and relatively small system sizes in that work raise ques- 
tions concerning its validity). If this scenario is correct, 
it suggests that the vortex glass transition line we find 
becomes only a crossover phenomenon in a real physi- 
cal superconductor with finite A. This crossover might 
still be observed in cur rent- volt age characteristics that 
seem to obey critical scaling, as at a second order glass 
transition, only to have the scaling break down at suf- 
ficiently small currents (which probe increasingly large 
length scales and thus are sensitive to finite A). The 
failure to conclusively demonstrate critical scaling in ex- 
perimental data that has been noted by Strachan et alpl 
may perhaps be due to such an effect. In contrast to 
this situation with respect to the vortex glass, we note 
that magnetic field fluctuations should not qualitatively 
effect the vortex lattice melting transition at low disor- 
der strength, since melting is believed to be mediated by 
short wavelength shear fluctuations and so should be rel- 
atively insensitive to behavior on the long length scale 
A. 

Looking at the vortex glass T g (p) and vortex lattice 
melting T m (p) transition lines in the T — p plane, shown 
in Fig. [I] it appears that these lines are of distinctly 
different origin, rather than one being a continuation of 
the other. However we are unable to directly investigate 
this issue as we are unable to sufficiently equilibrate our 
model system at the low temperatures where these two 
transition lines would appear to meet. 



thick as the size considered in section III , we find that as p 
increases above ~ 0.12 the sharp first order melting tran- 
sition that we observed in the thinner system now smears 
out over a finite temperature interval, T m 2 < T < T m i. 
In Fig. [7fi we show a plot of the ordering parameter 
AS /So vs simulation time for a particular realization 
of the quenched random disorder at disorder strength 
p = 0.18 and temperature T m2 < T = 0.18 < T ml . 
In contrast to what we found in Fig. [3] for the thinner 
system at weaker disorder, we now see that AS /So ap- 
pears to fluctuate about four different discrete values. 
In Fig. [TJd we plot the histogram P(AS) vs AS /So re- 
sulting from the data of Fig. [7^. We see clearly four 
distinct peaks which we identify as representing the lat- 
tice state (AS /So ~ 0.45), the liquid state (AS/S ~ 0), 
and what we will denote as the "mixed 1" and "mixed 
2" states (AS/S ~ ±0.2). The locations of the min- 
ima between these peaks, shown as the dashed lines in 
Fig. [T^i, we use as a simple criteria for sorting each of the 
microscopic configurations into one of these four different 
states. Having so sorted the microscopic states we can 
then compute the histograms of energy E and disorder 
conjugate variable Q for each of the four states. These 
we show in Figs. [TJ^ and d. We see that there is a discon- 
tinuous decrease in E, and a discontinuous increase in Q, 
as the system transitions from the liquid, to the mixed 
state, to the lattice. The distributions of E and Q are 
comparable for both mixed 1 and mixed 2 states. 

In Fig. [8] we show log intensity plots of the structure 
function, lnS'(k^), averaged separately over the configu- 
rations belonging to each of the four states. Figs. [8^l— d 
correspond to the lattice, liquid, mixed 1 and mixed 2 
states, respectively. We see that the lattice state, as ex- 
pected, has sharp Bragg peaks of relatively high intensity 
about the reciprocal lattice vectors Ki. The liquid state 
has only broad diffuse peaks of equal low intensity at 
both Ki and K 2 (see Fig.[2]3 for the definition of Ki and 
K 2 ). The mixed 1 and mixed 2 states have sharp peaks 
of intermediate intensity at Ki and K2 respectively. 

To investigate the nature of the mixed states, we 
plot in Fig. [9k the profile of S(k x ,k y ) vs k x for fixed 
k y = 4 (2ir/£) passing through the wavevectors Ki and 
K2. We see that the peak at Ki for the mixed 1 state of 
Fig. [8)3 is as sharp as, but smaller in amplitude than, the 
Bragg peak of the vortex lattice state of Fig. [8^1, while 
the remainder of this profile (in particular the small dif- 
fuse peak at K2) resembles that of the liquid state of 
Fig. |8b. A similar result holds for the mixed 2 state of 
Fig.|8p, expect that the behaviors at Ki and K2 are now 
interchanged. In Fig.|9]3 we plot the values of S(Ki) and 
5(K2), computed now for individual xy planes (instead of 
averaged over all such planes as in Eq. ([9])), vs the plane 
height z. We see that in the mixed 1,2 states roughly 



9 




0.08 

b) 

0.06 h q uid 

5? : 
< 0.04 

0.02 1 



# sweeps (x 10 6 ) 



T 



mixed 2 



0.00^ 



p = 0AS 7=0.18 
mixed 1 

lattice 





0.4 0.6 



0.12 r 



0.8 



1.0 



d) mixed 
0.08 liquid^ 

d 

oTO.04 
0.00 




-1.40 -1.39 -1.38 -1.37 -1.36 
E 



-0.40 



-0.30 -0.20 

Q 



a) ■ 








[ 




■ 






■ 



•k x 



c) . 






m. 






■ 




■ 






■ 










i 


■ 





b) 

■ 














FIG. 8: Logarithmic intensity plot of vortex structure func- 
tion lnS'(k^) averaged over the microscopic configurations of 
the (a) lattice, (b) liquid, (c) mixed 1 and (d) mixed 2 states, 
corresponding to the data of Fig. [7] 



FIG. 7: (Color online) a) Plot of AS = 5(Ki) - 5(K 2 ), nor- 
malized by So = S(k± = 0), vs Monte Carlo simulation time, 
at T = 0.18 and disorder strength p = 0.18, for a particu- 
lar realization of the quenched randomness in a system of size 
20 x 20 x 12. The unit of time on the plot represents 10 6 Monte 
Carlo sweeps through the system, and each data point repre- 
sents an average over 65536 consecutive sweeps. The system 
shows discrete transitions between states with four different 
average values of AS; these states are separated by the hor- 
izontal dashed lines and labeled "lattice", "liquid", "mixed 
1" and "mixed 2". b) Histogram P(AS) vs AS/ So for the 
data of part a, showing four peaks corresponding to the four 
different states. c,d) Histograms of energy E and disorder 
conjugate variable Q for each of the four states. 



half the system is ordered into a lattice, while half ap- 
pears disordered. Comparing the plots for mixed 1 and 
mixed 2, we see that the location of the more strongly 
ordered planes remains roughly the same, independent of 
the two possible orientations of the vortex lattice that is 
ordering in those planes. We therefore conclude that the 
relative order or disorder of planes arises from the specific 
distribution of random bond strengths {Ji^} that exists 
in the particular realization of the quenched randomness, 
rather than being an effect due to thermal fluctuations 
or being stuck out of equilibrium. 

To further understand the nature of the mixed state, 



in Fig. 10 we show intensity plots of the real space av- 
erage vorticity (n z (x,y, z)) in the xy plane for several 
representative layers at heights z. We show plots for the 
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FIG. 9: (Color online) a) Plot of vortex structure function 
S(k x ,k y )/So vs k x , for fixed k y = 4(2tt/L) passing through 
the reciprocal lattice vectors Ki and K2, for the lattice, liq- 
uid, mixed 1 and mixed 2 states, whose intensity plots are 
shown in Figs. [8^t— d. b)Plot of S(Ki) (solid symbols) and 
S(K.2) (open symbols), computed for the individual layer at 
height z, vs z. Results are shown for the lattice, liquid, mixed 
1 and mixed 2 states of Fig. [8] 
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FIG. 10: Intensity plots of real space average vorticity 
(n z (x,y, z)) in the xy plane at four representative layers, 
z = 0, 4, 7, and 10, for the lattice, mixed 1, and liquid states, 
corresponding to the data of Fig. [8] 



lattice, mixed 1, and liquid states, for the same run that 
produced the results shown in Figs. [7|-|9j A bright white 
square represents a pinned vortex that stays on that site 
throughout the course of the simulation run; a dark black 
square represents a site on which no vortex sits; squares of 
intermediate shades of gray represent sites that vortices 
hop into and out of during the course of the simulation. 
Hence a region of periodic white squares indicates a local 
region that has ordered into a pinned vortex lattice; a 
region of gray squares with no apparent structure to the 
relative shadings indicates a local region that is disor- 
dered into a vortex liquid. We see from Fig. [10] that even 
the state which is globally characterized as the lattice 
contains localized regions which are disordered; these are 
presumably regions where the vortex pinning, due to the 
local bond disorder {J^}, is stronger than average. In 
one layer, z = 7, we see large domains of each of the two 
different possible vortex lattice orientations (i.e. a do- 
main with ordering wavevector K2, as well as a domain 
with the dominant ordering wavevector Ki). In the state 
which is globally characterized as a vortex liquid, we still 
see individually pinned vortices. In the mixed state we 
see layers that are mostly lattice (z = 10), layers that are 
mostly liquid (z = 4), and layers that consist of coexist- 
ing domains of lattice and liquid (z = 0, 7). 



The above results suggest that the mixed state is one 
in which different regions of the system order from liq- 
uid to lattice at slightly different temperatures, due to 
fluctuations in local bond strength disorder. Depending 
on the particular realization of disorder, this can hap- 
pen in ways that are more complex than the particu- 
lar case illustrated in Figs, ff]— [10} In some other real- 
izations we have found the histogram P(AS) does not 
show such clearly separated peaks as in Fig. [7}}, and 
there may be suggestions of more than three peaks. In 
some realizations we find a mixed state in which some 
layers of the system have a lattice ordering specified by 
the wavevector Ki, while other layers of the system are 
ordered with the opposite lattice orientation specified 
by the wavevector K2. For such a case our ordering 
parameter AS = 5(Ki) — S^K^) can be very small, 
which by our previous criterion would misidentify such 
a state as a liquid. To avoid such misidentifications, we 
adopt instead a more general approach, looking at the 
two dimensional histogram P(S f (Ki), S^K^)). In such a 
histogram, the liquid state is represented by a peak at 
5(Ki) = S(K2) = 0, a uniformly ordered lattice state is 
represented by a peak at either 5(Ki) large and S^K^) 
small, or vice versa, and a mixed state is represented by 
a peak elsewhere in the S(Ki) x aS^K^) plane. In Fig.[TT| 
we show such a two dimensional histogram for a differ- 
ent particular realization of the randomness at disorder 
strength p = 0.14 and temperature T = 0.208. We see 
the lattice state with wavevector Ki coexisting with the 
liquid state, coexisting with a mixed state having roughly 
equal ordering in Ki and K2. Based on the location of 
the peaks in this two dimensional P(5'(Ki), S f (K 2 )) his- 
togram we adopt the following criteria for sorting mi- 
croscopic configurations into liquid, lattice and mixed 
states: if 5(Ki)/5 < 0.08 and S(K 2 )/S < 0.08 we 
assume the configuration is a liquid; if 5(Ki)/5o > 0.3 
and 5(K2)/So < 0-1 or vice versa, we assume the con- 
figuration is a uniform lattice; any other configuration is 
taken as belonging to the mixed state. 

As temperature is varied, the locations of the peaks 
in P(S f (Ki), S(K2)) vary only slightly, however their 
weights (total number of configurations in each peak) are 
found to vary in the following manner. As T decreases, 
the weight of the liquid state decreases while the weights 
of the other states increase. As T is decreased further, 
the liquid state disappears, the weight for the mixed state 
decreases, while the weight for the lattice state increases. 
Finally, at low enough T, the mixed state disappears and 
only the lattice state remains. We define T m i as the 
temperature where the liquid state and mixed state are 
roughly equal in weight. We define T m2 as the tempera- 
ture where the lattice state and mixed state are roughly 
equal in weight. Our determinations of T m i and T m 2 in 
this manner are made from rough eyeball estimates, as 
we find that the errors involved in such eyeball estimates 
are considerably smaller than the the sample to sample 
variation in these temperatures between different inde- 
pendent realizations of the quenched random disorder. 
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FIG. 11: Intensity plot of two dimensional histogram 
P(5 f (Ki,K 2 )) for a particular realization of the random dis- 
order at p = 0.14, T = 0.208. 



These values of T mi and T m 2, averaged over 8 indepen- 
dent realizations of the quenched randomness, determine 
the phase boundaries shown in Fig. [I] Choosing a tem- 
perature T m 2 < T < T m i where both liquid and lattice 
coexist, we also compute the averages of E and Q in the 
liquid and lattice states separately. We then compute AE 
and AQ as in Eq. (15), representing the total change in 



these quantities as the system orders from liquid to lat- 
tice. For a few random realizations, where our simulation 
temperatures did not include a value where liquid and 
lattice coexisted, we used standard extrapolation meth- 
ods to extrapolate the liquid (lattice) histogram from a 
slightly higher (lower) temperature to a common tem- 
perature T m2 < T < T m i, and then computed AE and 
AQ from these extrapolated histograms. Our results, 
averaged over 8 independent realization of the quenched 
random disorder, are shown in Fig. [5|i (solid symbols) 
where we see that they appear to be a consistent con- 
tinuation of the results found in the smaller size system 
(20 x 20 x 6) at weaker disorder strengths p, except that, 
as noted previously in Table [l[ [AQ] is somewhat larger 
in the larger system size. 

The above results suggest that as the system size and 
the disorder strength increase, the transition from vor- 
tex liquid to lattice remains locally first-order-like, as 
was found in Section III, however different regions of 



the system may order at slightly different temperatures. 
Such a conclusion is in perfect agreement with experi- 
ment results from magneto-optical imaging of vortex lat- 
tice melting in BSCCO, 27 where it was observed that the 
spread in local melting temperatures over the area of the 
sample increased as the average melting temperature de- 
creased (i.e. as the average pinning strength increased). 
A similar observation of coexisting ordered and disor- 
dered vortex regions over a finite temperature interval at 
the vortex lattice melting transition was reported in more 



recent experimental studied of the peak effect region of 
NbSe 2 . 

The vortex lattice melting transition thus becomes 
smeared out over a range of temperatures when viewed 
on the global scale. Such a scenario was predicted long 
ago by Imry and Wort is who considered the effect of 
quenched randomness on a first order phase transition in 
terms of a competition between the decrease of the free 
energy due to local ordering of domains of finite size in 
regions of lower than average quenched disorder, vs the 
increase of free energy due the surface tension of the re- 
sulting domain walls enclosing the ordered domains. We 
have attempted a quantitative test of the Imry-Wortis 
scenario as applied to our model. Our efforts are de- 
scribed in Appendix A. While they are suggestive, show- 
ing the right trends as disorder strength p increases, our 
results are not conclusive. Even if the Imry-Wortis sce- 
nario applies to our mixed state, two possibilities still 
exist: (i) it may be that on the global scale in the ther- 
modynamic limit, there remains a sharp first order phase 
transition with reduced but still discontinuous jumps AE 
and AQ at a single well defined T m , or (ii) it may be that 
there will be multitude of local transitions spread out 
smoothly over a range of spatially varying local melt- 
ing temperatures, thereby converting the transition on 
the global long length scale to a continuous second order 
transition. To discriminate between these two possibili- 
ties would require looking at much larger systems sizes 
than we are currently able to equilibrate. 



VI. DISCUSSION 

To summarize, our results for the vortex phase dia- 
gram of the uniformly frustrated 3D XY model with dis- 
ordered couplings show many qualitative features in good 
agreement with experiments on strong type-II supercon- 
ductors. We find a sharp local first order vortex lattice 
melting transition at low pinning strengths p. The melt- 
ing temperature decreases as p increases, and appears 
to be steadily decreasing towards zero as a critical pin- 
ning strength p c is approached. We find that this melting 
transition remains first order, with no evidence for a crit- 
ical end point or other multicritical points, down to the 
lowest temperatures to which we can equilibrate. At high 
pinning strengths p > p c we find in our model a sharp 
second order vortex glass transition. This glass transition 
may evolve into a non-signular crossover phenomenon if 
magnetic field fluctuations, due to a finite magnetic pen- 
etration length A, were incorporated into the model. A 
completely new feature of our simulations is the observa- 
tion that at intermediate disorder strengths p < p c the 
vortex lattice melting transition is smeared out over a 
temperature interval of finite width, corresponding to co- 
existing regions of ordered and disorde red v ortex states, 
as has been seen in recent experiments! 7 -^ 

It is interesting to compare our results to those of 
Nonomura and Hu 45 (NH), who studied a very similar 3D 
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XY model but with much more weakly coupled planes, 
J z /J± = 1/400, a more dilute vortex density / = 1/25, 
and with random couplings that model a dilute set of 
localized strong pins rather than the amorphous random 
couplings we have used here. In addition to vortex lattice 
melting and vortex glass transitions, NH reported the ex- 
istence of a vortex slush phase, lying between the vortex 
liquid and the vortex glass, and separated from the vortex 
liquid by a sharp first order transition. In earlier worlP^ 
we have repeated simulations of Nil's model, using the 
exact same parameters as NH. We find that their vortex 
slush phase shares some similarities with our interme- 
diate region discussed above, in that both are regions in 
which the vortex lines have only partially ordered. There 
are however some important differences. 

(i) Our intermediate region lies between the vortex liq- 
uid and vortex lattice phases and so may be thought of 
as a broadening of the melting transition; when we cool 
at fixed p though the intermediate region, our system 
(except in rare cases when we fail to equilibrate) always 
orders into a clear vortex lattice. In NH's model, the 
vortex slush lies above the vortex lattice phase (similar 
to what was reported in some experiment J 9 * 26 ^ ) ; cooling 
through NH's vortex slush, one enters the vortex glass 
and not the vortex lattice, (ii) We showecP^J that con- 
siderable hysteresis existed in the region of NH's vortex 
slush phase. In our present simulations there is no hys- 
teresis: at fixed simulation parameters in our interme- 
diate region the system is repeatedly hopping into and 
out of the vortex lattice, liquid, and mixed states (see 
Fig. [7^i) and our system is thus fully equilibrated, (iii) 
We argued 58 that in NH's model, most planes of their 
vortex slush contained an ordered vortex lattice, however 
the orientations of the vortex lattice varied with height z\ 
we argued that these mismatched vortex lattice orienta- 
tions would be unfavorable in the thermodynamic limit 
and that NH's vortex slush was most likely a finite size 
effect to be replaced by an ordered vortex lattice as sys- 
tem size increased. In our intermediate region, however, 
we see coexisting planes of mostly vortex lattice, planes 
of mostly vortex liquid, as well as planes with large do- 
mains of both lattice and liquid. The scaling argument 
we used against NH's vortex slush thus does not apply. 
It is of course possible that upon increasing system size, 
NH's vortex slush will similarly develop coexisting or- 
dered and disordered domains within individual planes 
and so remain as a stable phase. 

We conclude therefore that our intermediate region is 
distinctly different from the vortex slush of NH, and is 
perhaps more similar to the multidomain glass state that 
has been proposed by MenonP Unfortuntely, we have 
not been able to equilibrate our system in the interesting 
region where the vortex lattice melting transition meets 
the vortex glass transition. 

Another set of interesting simulations has been carried 
out recently by Dasgupta and Vali d 5 1 | 52 J using a density 
functional approach applied to interacting pancake vor- 
tices in a 3D layered system. They consider both the case 



of dense amorphous pins J 5 ^ such as we consider here, and 
the case of dilute well localized pinsj^ closer to the model 
of NH, mapping out the phase diagram as a function of 
temperature and pinning strength. Their approach, be- 
ing essentially an equilibrium mean field method in the 
presence of quenched randomness, is suited to locating 
first order phase transitions, as in vortex lattice melting 
or the proposed vortex slush, rather than continuous sec- 
ond order transitions, as one expects for a vortex glass. 
For the dense amorphous pinning, they find only a sin- 
gle unified vortex lattice melting transition, qualitatively 
similar in shape to what we find in the present work. 
Their vortex liquid state shows no significant local order- 
ing on length scales larger than the average vortex spac- 
ing and the average vortex density varies smoothly as one 
goes from the vortex liquid at weak pinning to the vor- 
tex liquid at strong pinning. They find no vortex glass, 
no vortex slush, and no multidomain glass of polycrys- 
taline domains. For the case of dilute well localized pins, 
they find again a vortex lattice melting transition with a 
similar shape as before (though quantitatively at a very 
different location, comparing amorphous to dilute pins 
with equal second moments of the random pinning po- 
tential). Their vortex liquid state, however, now shows a 
clear polycrystaline structure with noticeable short range 
translational order extending on lengths larger than the 
average vortex spacing. However they again find no first 
order transition within their vortex liquid phase, such 
as might define a region of vortex slush or multidomain 
glass as distinct from the vortex liquid. Using a perco- 
lation criterion to define a crossover to glassy behavior 
within the vortex liquid phase, they find a line that is 
qualitatively similar in location to our vortex glass tran- 
sition, as shown in Fig. [I] 
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Appendix A 

We give a simplified summary of the Imry-Wortis sce- 
nario, as applied to our model, as follows. Consider T m (p) 
the nominal melting temperature of a system with disor- 
der strength p. Let 

A/(p, T) EE /latticed, T) ~ /liquid T) (19) 

be the difference in free energy density between the vor- 
tex lattice and liquid states. Consider now a volume v 
in which, due to the random distribution of pins, the 
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FIG. 12: 

effective disorder strength p e ff is either greater than, or 
less than, the average p, p e ^ = p ± Ap (Ap > 0). Since 
dT m /dp < 0, if Peff > the domain would lower its bulk 
free energy by disordering for some range of temperatures 
AT below T m (p). Similarly, if p e ^ < p, the domain would 
lower its bulk free energy by ordering for some range of 
AT above T m . This is sketched in Fig. [12] below. For 
the system siting at T m (p), the decrease in the bulk free 
energy of such a domain would be, 

AT bulk = ±vAf(p±Ap,T m (p)) 
« ^A/(p,T m (p)) ^ 
dp 

= vAQ(p,T m (p))Ap . (20) 

where the last line follows from our definition of AQ in 
Eq. (15). 



On the other hand, the formation such a domain would 
lead to an increase in surface free energy due to the re- 
sulting domain wall, 



AT surface - 2aJ 2 ± + Aa ± £ ± £ z 



(21) 



where we take account of the anisotropy of the system by 
denoting £ z and £± as the lengths of the domain parallel 
and perpendicular to the applied magnetic field, with a z 
and ct_l the surface tension between vortex lattice and 
liquid for surfaces with normal in these respective direc- 
tions. The domain will be unstable to flipping only if the 
total change in free energy is negative, 



AT = -ATbuik + AT sur f ace < . 



(22) 



For simplicity, we will assume that the domains which 
form are such that the surface tension is equally dis- 
tributed over all surfaces, so that, 



so £ z 



(23) 



and so the volume of the domain is v ~ £\£ z 



(o~ z /cr±)£j_. The instability condition Eq. (22) then be- 
comes, 



AF(£ ± ) 



0~± 



£ 6 ± AQAp + 6o- z r ± < . 



(24) 



Next, we can write that the typical variation in disorder 
Ap for a domain of size £\£ z can be written as, 

AT m AT m A£ 
* P \dT m /dp\ T m AQ ' 



where AT m is the variation in melting temperatures for 
domains of size £\£ z sampled from a system with av- 
erage disorder p, and the last equality follows from the 
Clausius-Clapeyron relation, Eq. (18). If the system is 
self averaging, as is suggested by our results in Table [T| 
then we expect 



AT m 



a 



(26) 



where a is some constant. Substituting Eqs. (26) and 



(25) into (24) then gives 



AF(£ ± ) = 



r o T z ~ aAE 3/2 



6a z £\ < . 



(27) 



We sketch AF(£ ± ) in Fig. [13 
AF(£ ±0 ) 



Defining £±q by 



0, we see that the system can be unstable 
only to the flipping of domains with transverse length 
£± < £±o- From Eq. (f27|) we have, 



1 



t_L0 



o~ z o~± 



aAE 



(28) 



However, the above arguments only apply to domains 
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which are well defined as such, i.e. they are at least as big 
as the correlation length that defines the minimum size 
of a domain. Hence the instability condition for domain 
flipping becomes, 



(29) 



and so, in particular, the system is stable against the 
flipping of local domains if £_l > £_\_o- As the disorder 
is increased, one in general expects £j_q to increase. So 
if at low disorder the system is stable, £j_ > £±o, as the 
disorder increases one will eventually reach the condition 
£x = £±o an d the system will first become unstable to 
domains on the size of the correlation length As the 
disorder increases further, larger domains of size £±, with 
> £± > will go unstable. 
To test the Imry-Wortis scenario for our vortex line 
system, we therefore wish to compute the length £±q of 
Eq. (28). We have already computed T m and AT, as 
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shown in Figs. [T] and [5] respectively. We use our results 
from the 20 x 20 x 6 size system, averaging AE/T m over 
the different realizations of randomness, and computing 
a = AT m ^/V from the observed spread in melting tem- 
peratures, such as shown in Table [T] for the specific case 
of p = 0.12. It remains, therefore, to compute the surface 
tensions <j z and a±. 

To compute the surface tension we use a method based 
on the approach of Potvin and RebbiP^ We take a given 
realization of the randomness for which we have previ- 
ously determined the melting temperature T m . We then 
take an exact copy of this system and join it to the orig- 
inal along the surface whose surface tension we seek to 
compute. On one side, denoted as "side 1", we use cou- 
plings J_li = Jj_(l + Si) and J z \ = J z (l + Si) while 
on the other side, denoted as "side 2" , we use couplings 
Jj_2 = + S2) and J Z 2 = J z {\ + S2). In this way we 
expect that exactly at T m (as determined in the origi- 
nal system with Si^ = 0) if Si^ > 0, that side will be 
ordered, while if Si^ < 0, that side will be disordered. 
Choosing Si = —S and S2 = +S will thus create an inter- 
face between ordered and disordered halves of the total 
system. Consider now a trajectory in the (Si ^2) plane, 
as shown in Fig. 14 In this figure, point A is a totally 




above equation we use, 

dF(S u S 2 ) 



F A 



+<5 



dSi - 



dSi 



-s 



S 1 



(32) 

where Ei(Si, S2) is the total energy of side 1 at the speci- 
fied couplings. A similar expression can be derived for 
Fb — Fc- Simulating at points along the trajectory 
A — > B — > C we then integrate the energies Ei and E2 
to compute the surface tension, 



E 
A 



1 

AA 



+5 



dSi 



E 1 (6 1 ,S 2 ) 

1 + 8! 



+S ,, E 2 (8!,5 2 ) 

-S 1 + 02 



(33) 

where A is the total area of one interface. We implement 
this procedure on a 20 x 20 x 6 system doubled in the z 
direction (to make a 20 x 20 x 12 system) so as to compute 
cr Zl and doubled in the x direction (to make a 40 x 20 x 6 
system) so as to compute cr±. We use a value S — 0.1 in 
order to get reasonable results. Our results are averaged 
over 8 independent realizations of the random disorder 
(only 7 for p = 0.18). We plot our results in Fig. 15 As 
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disordered system, point C is a totally ordered system, 
and point B has side 1 ordered and side 2 disordered. 
The total free energy of the system at point B can be 
written as 

F B = ^F A + iF c + 2E (30) 

where Fa and Fc are the total free energies at points 
A and C, and E is the total surface free energy of one 
interface between the ordered and disordered halves. The 
factor 2E appears since our periodic boundary conditions 
necessarily creates two interfaces equally spaced by half 
the length of the total system. From this we have, 

4E = [F B - F A ] + [F B - F c ] . (31) 

The surface tension between coexisting disordered and 
ordered phases at the same transition temperature T m 
is then obtained from E, taking in principle the limit of 
S — > 0. To evaluate the free energy differences in the 



FIG. 15: Surface tensions between vortex lattice and vor- 
tex liquid states, for surfaces with normal parallel (a z ) and 
perpendicular (a±) to the applied magnetic field, vs average 
disorder strength p. 

expected, a z and a±_ decrease as the disorder strength p 
increases. For our parameters of anisotropy and vortex 
line density we find <j z ~ <J±/3. 

We summarize the pieces of our calculation of £±0 in 
Table [TTJ The values for a z and a_\_ are obtained as de- 
scribed above. Values for [AE/T m ] are obtained averag- 
ing over careful equilibrations of 20, 8, and 7 different 
realizations of the random disorder for p = 0.12, 0.16 
and 0.18 respectively, for a 20 x 20 x 6 system. Because 
the spread in melting temperatures AT m is the quantity 
that is most sensitive to the fact that we sample only 
over a rather small number of random realizations, for 
p = 0.16 (0.18) we have tried to do better than the 8 (7) 
realizations we have carefully equilibrated by computing 
AT m from 16 random realizations where we determine 
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T m from shorter runs and more qualitative methods. We of Fig. [T] shows that the mixed state is already observed 
then use a = AT m V20 2 x 6. for disorder strengths as low as p — 0.14. 



TABLE II: Values that enter our calculation of 



from 



Eq. (28) 



p 


AT m 


[AE/Tm] 


cr z 


o~± 


£±o 


0.12 


0.0037 


0.124 


0.0051 


0.0153 


0.18 


0.16 


0.0060 


0.074 


0.0034 


0.0091 


0.42 


0.18 


0.0090 


0.076 


0.0025 


0.0051 


2.41 



Our results show £±o to be an increasing function of 
disorder strength p : as expected. However, for the sys- 
tem to become unstable to the flipping of domains it 
is necessary that £±o > For our vortex density of 
/ = 1/5 we estimate that £_l at T m is at least as large 
as the average vortex spacing, a v = l/y/E ~ 2.2. This 
seems consistent with the real space images of Fig. 10 
were we see ordered regions of at least this size in the 
liquid, and disordered regions of at least this size in the 
lattice. Thus we have £j_q > £_l only for our strongest 
disorder strength p = 0.18, where our results are per- 
haps the least accurate. In contrast, our phase diagram 



It should be noted that the above analysis is based only 
on typical "root mean square" behaviors. Correlations 
between bulk and surface free energies of domains may 
enhance the effect over what we have estimated above. 
For example, a domain of lattice may flip to the liq- 
uid state in a region where the vortex pinning is locally 



stronger than on average; but in such a region, Fig. 15 
shows that the surface tension is lower than average, thus 
reducing the energy cost of such a flip from that consid- 
ered in our arguments above. Domains may also flip in 
regions of the system where the value of the local disor- 
der strength lies further out in the tails of the disorder 
strength distribution, rather than near the root mean 
square value. This might explain why we see our inter- 
mediate mixed states more easily when we increase the 
system size, thus affording a wider sampling of the disor- 
der strength distribution within any given single sample. 
The results of our Imry-Wortis analysis thus show the 
right trends for explaining our intermediate mixed states, 
however in the absence of a clear quantitative agreement, 
£±o ~ (i, we must regard our results as still inconclusive. 
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